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Abstract: One of the fundamental aspects of statistical be- 
haviour in many-body systems is exponential divergence of neigh- 
bouring orbits, which is often discussed in terms of Liapounov 
exponents. Here we study this topic for the classical gravi- 
tational iV-body problem. The application we have in mind 
is to old stellar systems such as globular star clusters, where 
A'" ~ 10^, and so we concentrate on spherical, centrally concen- 
trated systems with total energy E < 0. Hitherto no connection 
has been made between the time scale for divergence (denoted 
here by tg) and the time scale on which the energies of the par- 
ticles evolve because of two-body encounters (i.e. the two-body 
relaxation time scale, tr), even though both may be calculated 
by similar considerations. 

In this paper we give a simplified model showing that di- 
vergence in phase space is initially roughly exponential, on a 
timescale proportional to the crossing time (defined as a mean 
time for a star to cross from one side of the system to another). 
In this phase te <^ if N is not too small (i.e. N » 30). Af- 
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ter several e-folding times, the model shows that the divergence 
slows down. Thereafter the divergence (measured by the ener- 
gies of the bodies) varies with time as t^/^, on a timescale nearly 
proportional to the familiar two-body relaxation timescale, i.e. 
te ~ tr in this phase. These conclusions are illustrated by nu- 
merical results. 

Keywords: Gravity, Few-body systems. Relaxation processes. 
Particle orbits 

1. Introduction 

The classical gravitational A^-body problem is defined by the 
equations 

r, = -G E m,^^ (1) 
j + i 

where is the three-dimensional position vector of the i^^ star, 
rrii is its mass, and G is the universal constant of gravitation. 
We consider apphcations in which the total energy, E, in the 
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barycentric frame is negative and the total angular momentum is 
negligible. Starting from a rather broad set of initial conditions, 
such solutions settle down into a roughly spherical distribution 
of bodies in approximate "dynamic equilibrium" (Fig. 1), i.e. 
the spatial distribution is nearly time-independent on the time 
scale of the orbital motions of the particles. 

Early numerical integrations'®' with < 32 showed that 
a small change in initial conditions led to a roughly exponen- 
tial divergence of solutions (measured in 6A'"-dimensional phase 
space), even though the spatial distribution of the bodies in 
the two solutions might be indistinguishable within statistical 
fluctuations. The timescale of divergence, tg, was of order the 
crossing time, tcr, defined in a certain conventional way as the 
time for a body with a typical speed to move a distance compa- 
rable to the size, R, of the spatial distribution of the particles'^!. 
Thus 

R , , 

where V is the root mean square speed of the particles. Later 

3 



work'^'^'^l extended numerical results to larger N, and Good- 
man et al'"^l devised theoretical models confirming that te/tcr is 
virtually independent of N. 

One particular statistical specification of the initial condi- 
tions which has been studied is the Plummer model, which is 
often used in stellar dynamics for the study of relaxation and 
related processes. It is the stellar dynamical analogue of the 
n — 5 polytrope. For this model it has been found^ that 

" ~ In (0.73 In AT)- 

The functional form is suggested by a theoretical model '^1 , and 
the coefficients are not thought to depend sensitively on the 
initial conditions. Therefore for large star clusters generally we 
have 

te ~ O.Odtcr- (3) 
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Fig.l. Spatial distribution of bodies in a typical simulation. On 
left is a snapshot, and on the right is the numerically generated 
space density as a function of radius. 

The theoretical models of Goodman et al'^l dealt with the 
linear divergence of neighbouring solutions, when the separation 
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in position of the i body satisfies the variational equation 
■• K f Sri- Srj {5ri - 5rj) . (r^ - r^) \ 

(4) 

For practical purposes, however (e.g. for understanding the 
growth of errors in a numerical integration) the resulting roughly 
exponential growth quickly leads to separations so large that the 
linear approximation fails. In this contribution we develop the 
simplest model of divergence to account for the later, nonlinear 
growth of the separation between neighbouring solutions. We 
shall see that the time dependence changes from roughly expo- 
nential to roughly power-law, and that the timescale changes 
from roughly the crossing time to nearly the two-body relax- 
ation timescale, tr- This is the timescale on which the energies 
of the individual bodies vary significantly. Standard theoryt^'^1 
shows that 
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for systems of the general kind considered here. 

2. A model of divergence 

2.1 Linear growth of errors 

In this section we introduce a toy model for the divergence 
of neighbouring orbits. Though it gives much insight into the 
physics of the problem, many details are omitted. In the first 
instance we apply it to the linear regime in which the approx- 
imate eqs.(4) are valid. In this regime more elaborate models 
have been constructed by Goodman et al'^l . 

We make the following assumptions. As in the theory of two- 
body relaxation'^'^] we assume that the trajectory of a particle 
is nearly rectilinear, except for occasional two-body encounters 
(Fig. 2). We suppose that the important encounters are in the 
small-angle scattering regime, such that p 3> Gmjv^ where p 
is the impact parameter and v is the relative velocity of the 
two particles. In computing the effect of one encounter, we 
suppose we can treat the scatterer as fixed. We also suppose 
that successive encounters can be treated as if motion takes 
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place on one plane, and that the difference between two orbits 
is measured by the difference in the impact parameter, 5p. We 
assume that all particles have the same mass m. Finally, we 
suppose that the system is in virial equilibrium (see Binney & 
Tremaine^), which implies that 



GmN , , 

V __. (6) 



Here the symbol ~ means "is of order" , i.e. that the relation 
is approximate, and any numerical coefficient is ignored. Thus 
V r^V ^ for example. 




Fig. 2. Two successive encounters. 



In the small-angle scattering regime the maximum accelera- 

CjTft 

tion of the moving particle is of order — — and the duration of 
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the encounter is of order — . Thus the change in velocity is of 

order ^IIL^ an^l so the angular deflection is of order — - (Fig. 2). 
pV pV^ 

After the scattered body has travelled a further distance D to 

■ . . n . . r 1 GmD 
its next encounter, its spatial deflection is of order 



pV^ ' 

Now suppose the body had approached the first encounter 

on a parallel path at a slightly different impact parameter p+5p. 

Then, at the time of the second encounter, its position would 

GmD 

have been displaced by a distance of order 6p H — — The 

first term is the displacement that would have occurred even 
in the absence of the first encounter. The second occurs be- 
cause, if 5p > 0, the body has been deflected less by the first 
encounter. (The differential approximation used for this term 
is not valid unless \6p\ <^ p; this is the approximation which 
restricts the present theory to the linear regime in which eqs.(4) 
are valid.) The total displacement measures the change in im- 
pact parameter at the second encounter. Hence the variation in 

p is multiplied by a factor of order ^1 H TyY^ P^'^ encounter^. 

^In a fully three-dimensional treatment this becomes a matrix equation. 
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Now we consider the cumulative effect of several encounters 

within a restricted range of impact parameters around the value 

p, e.g. from p/2 to 2p, but ignoring other encounters. We 

start at some time t and consider the effect of encounters in 

a subsequent interval At, chosen sufficiently large that several 

encounters occur within this interval. The actual number of 

AtV 

such encounters is of order ^ , and so the variation in the 
orbit is givenQ by 

r / X A GmD\ D 

Also it is clear that p^Dn ~ 1, where n is the number of particles 
per unit volume, and so 

5r(t + At)~5r(t)(^l + -j^j . (7) 

It follows from the relation n ~ — and eq. (6) that 



Sr {t + At) ~ 5r (t) 1 + 



AtNp^ 



^We ignore two complications which tend to counteract each other: (i) 
the persistence of effects of early encounters, and (ii) partial cancellation 
of successive encounters by their vectorial character. 
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where we have used eq.(2). 

Encounters take place at a wide range of impact parameters 
p. Writing eq. (8) as 

In 6r {t + At) - In 5r (t) ^ -N^^ In (^1 + -^^ j (9) 

we see that those with p « RN'^^"^ are individuaUy very effec- 
tive but too rare to dominate, whereas those with p» RN~^^'^ 
lose out by being individually ineffective, despite being very nu- 
merous. Encounters at impact parameter p ~ RN'^/"^ are 
most effective cumulatively, and lead to exponential growth of 
the deviation 5r, on a timescale of order tcr- 

Another way of seeing this is to sum the right hand side 
of eq.(9) over all impact parameters p <j R. Since this term 
represents the effect of encounters with impact parameters in a 
range near some value p, the summation can be accomplished 
by multiplying by dp/p and integrating. The result is that 

In 5r {t + At) - In 5r (t) ~ — In A^. 

tf-r 

Except for the logarithmic factor, this is equivalent to the re- 
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suit obtained by ignoring all encounters except those near p ~ 

RN-^'\ 

Many factors have been omitted from this simple model, 
including the distribution of velocities and density, and the 
curved orbits of bodies between encounters. Nevertheless, the 
results of more detailed models and numerical simulations, al- 
ready quoted, confirm our basic result, except for a very weak 
A'"-dependence. 

2.2 Nonlinear growth of separation 

The above theory is valid as long as 5r « p, and here 
we may take for p the impact parameter for the most effective 
encounters, i.e. p ~ RN~^/'^. Suppose we are interested in 
growth of errors in an A^-body integration of eqs. (1), for a 
system which has been scaled so that ~ 1. Then we may 
have 5r(0) ~ 10~^^ for a double precision calculation, and so the 
linear approximation breaks down after about 30ie, i-e. between 
one and two ten by eq.(3). 

Thereafter we suppose that encounters with impact param- 
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eters p << Sr are ineffective. Then we may estimate the 
growth of the separation of neighbouring orbits by substitut- 
ing p ~ Sr{t) in eq. (9), which gives 

\nSr{t + At) -In Sr{t) r ^^in 1+ > 



We are in a regime where 5r{t) ^ RN~^/'^, and so we can ap- 
proximate 

luM< + A*)-lnM«)~^^. (10) 

Since the term on the right depends on t, we can no longer 
conclude that In 5r{t) increases linearly with t. To determine its 
time dependence we rewrite eq. (10) as a differential equation, 



I.e. 



dt "tcr N5r{t)^ 
Ignoring for the moment the distinction between "~" and , 

we obtain the solution 

where to is a constant, which may be interpreted as the time at 
which the growth of errors enters the nonlinear regime. 

13 



Well into the nonlinear regime we now see that Sr (t) ~ 

/ t \ 

R y jy^ J ■ III order to interpret this result we shall estimate 

the difference in binding energy, e, of the body between the two 

GNm 

neighbouring solutions. Now e ~ — - — , and we can estimate 

R 

G NmSr 

Se ~ — — — . (We could obtain a similar estimate from con- 

Se ( t 

sideration of the difference in velocity.) Hence — ' ' 



Now the two-body relaxation time, tr, may be estimated by 
eq.(5), and so — ~ ( ~ ) if we ignore a logarithmic depen- 
dence on N. 



3. Discussion 

Recall that we are considering two solutions of eq. (1) start- 
ing with shghtly different initial conditions. Suppose that we 
measure the separation of the two solutions by the separation 
in energy, de, of a typical body. What we have concluded is 
that, for at most a few crossing times, 5e{t) grows exponentially, 
with an e-folding time comparable with tcr itself. Thereafter 
Se{t) approaches a power law dependence, varying as t^^^, on a 
timescale of the relaxation time. 
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The standard theory of relaxation tells us how e (the energy 
of a given star) evolves on a single solution of the A^-body equa- 
tion. If we ignore variations of e inside an encounter, e performs 
a random walk on the timescale tr, and the change in e varies 
as t^/^. (We here ignore the role of "dynamical friction", which 
corresponds to the drift term in a Fokker-Planck description of 
the relaxation'^'^1.) 

Fig. 3 illustrates these points using data from numerical N- 

body integrations with = 256. Two systems were integrated 

simultaneously using identical initial conditions except for a 

small difference in one coordinate of one particle. The solid 

curve (a) shows the mean square difference in the energies of 

the N particles^. The corresponding initial conditions were also 

used for simultaneous integration of the variational equations, 

and the long-dashed curve (b) shows the corresponding mean 

square variation of energy. This grows nearly exponentially, but 
■^Similar results have been presented by Merritt^^l for motion in the 

gravitational field of N fixed bodies 
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is followed by (a) for only a limited time of order a crossing 
time. The short- dashed curve (c) shows the mean square differ- 
ence between the initial energy and the energy at time t, again 
averaged over the N particles. This is caused by two-body re- 
laxation. Evidently curve (a) departs from curve (b) around the 
point where the latter crosses curve (c) , and then nearly follows 
(c). In this way we see that the growth of errors, which is expo- 
nential only in the linear regime, is consistent with the theory 
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of two-body relaxation. 
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Fig. 3. Mean square energy difference in numerical integra- 
tions with N — 256, as a function of time. The meaning of the 
different curves is stated in the text. The results plotted are 
the mean of four independent runs. In the adopted units the 
crossing time is 

The exponential divergence slows down to a power- law growth 
because close encounters become increasingly ineffective. There 
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is a geometric way of looking at this. Krylov'^l showed that 
the divergence could be understood as the behaviour of neigh- 
bouring geodesies on a certain manifold. As two neighbouring 
geodesies deviate further, their deviation is influenced less and 
less by the fine geometrical structure of the manifold across 
which they are proceeding. 
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